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ABSTRACT 

In this paper, gravothermal oscillations are investigated in multi-component star 
clusters which have power law initial mass functions (IMF). For the power law IMFs, 
the minimum masses (m m i n ) were fixed and three different maximum stellar masses 
{m max ) were used along with different power-law exponents (a) ranging from to 
—2.35 (Salpeter). The critical number of stars at which gravothermal oscillations first 
appear with increasing N was found using the multi-component gas code SPEDI. The 
total mass (M to t) is seen to give an approximate stability condition for power law IMFs 
with fixed values of m max and m min independent of a. The value M tot /m max ~ 12000 
is shown to give an approximate stability condition which is also independent of m max , 
though the critical value is somewhat higher for the steepest IMF that was studied. 
For appropriately chosen cases, direct N-body runs were carried out in order to check 
the results obtained from SPEDI. Finally, evidence of the gravothermal nature of the 
oscillations found in the N-body runs is presented. 

Key words: globular clusters: general; methods: numerical; methods: n-body simu- 
lations. 



1 INTRODUCTION 

The condition for the onset of gravothermal oscillations is 
best understood for the case of one-component star clusters, 
clusters consisting of stars of equal mass. [Goodman (19871 



found that gravothermal oscillations first appear when the 
number of stars N is greater than 7000. This condition has 



also been confirmed with Fokker-Planck calculations ( Cohn 
|et al|1989| and by direct N-body simulations ( |Makino|1996[ ). 
However, the multi-component case is more complicated. 
This is due to the fact that the presence of several compo- 
nents introduces new dynamical processes into the system, 
and several additional parameters in addition to N. 

Even for the two-component case, which is the sim- 
plest kind of mass spectrum, the condition for the onset of 
gravothermal oscillations is not so simple. Two-component 
models can be subdivided into Spitzer stable and Spitzer 
unstable cases depending on whether or not the two compo- 
nents can achieve cquipartition of kinetic energy during core 
collapse ( |Spitzer||1987[ ). |Kim, Lee fc Goodman| \1998\ stud- 
ied a range of Spitzer stable two-component models. Their 
research supported the applicability of the Goodman stabil- 
ity parameter e (see Goodman 1993) as a stability criterion. 



Breen & Heggic (20121, whose research focused on the more 
general Spitzer unstable two-component case, indicated that 
the occurrence of gravothermal instability depends approx- 
imately on the number of stars in the heavier component. 
Breen & Heggie (20121 also found that the critical value 
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of e depended on the parameters of the mass function (e.g. 
stellar mass ratio). However, by using a slightly modified 
version of e, one with a modified definition of the half mass 
relaxation time, they found a nearly constant critical value. 

|Murphy et al| ( |1990| ) found that the post-collapse evolu- 
tion of multi-component models was stable to much higher 
values of N than in one-component models and that the 
value of N at which gravothermal oscillations appeared 
varied with different mass functions. They studied seven- 
component systems constructed to approximate evolved 
power law IMFs with masses ranging from 0.1 to 1.2Mq. 
The power law exponent that they considered ranged from 
—2 to —4.5. They found that gravothermal oscillations ap- 
peared when the to tal mass of the syste m (Mtot) was of or- 
der 8 x 10 4 Mq (see Murphy et al|l990| Figure 6) and that 
the critical value of Mtot increased with decreasing power 
law exponent. They suggested that the appearance of oscil- 
lations depends on the number of heavier stars. However, 
this leads to the issue that in a multi-component system it 
is not clear what the definition of a heavy star should be 
(this point is discussed in Section 2). 
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The main aim of the present paper is to provide a 
theoretical understanding of the onset of gravothermal os- 
cillations in multi-component systems. As this present pa- 



per follows on from the research of Breen & Heggie (20121 
it is worthwhile attempting to extend the concepts devel- 
oped in that paper to the multi-component case. Although 
two-component systems may be realistic approximations of 



Table 1. Critical value of N (N cr n) in units of 10 4 . The values 
of N cr i t in brackets were obtained using 5-component models, 
while all other values were obtained using 10-component models. 



The value of N e j for the case a 



-2.35 and extreme masses 
(3,0.1) could not be obtained with a 10-component model due to 
numerical difficulties. 



multi-component systems (Kim & Lee 19971, it is best to 
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multi-component systems as real globular clusters contain a 
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is the effect of varying the maximum stellar mass (m max ) on 
the onset of instability as this was not studied by |Murphy| 
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et al (19901. 
















The rest of this paper is structured as follows. In Section 
















2, we state the results concerning gravothermal oscillations 


Table 2. Critical value of M (M crit ) in 


units of 10 




in gaseous models. This section also contains subsections on 
















the Goodman stability parameter and a variant which used 
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a modified relaxation timescale. This is followed by Section 3 
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in which the results of N-body simulations are given. Finally, 
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Section 4 consists of the conclusion and discussion. 
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2.2 Interpretation of results 



Guided by the results of Murphy et al ( 1990 1, we first con- 



CRITICAL VALUE OF N 



2.1 Results of gaseous models 

In all cases, the initial conditions used were realisations of 



the Plummer model (Plummer 1911 Heggie fc Hut 20031 



The initial velocity dispersion of all components and the 
initial ratio of density of all components were equal at all 
locations. The initial conditions were constructed in order 
to approximate a continuous power law IMF with different 
exponents a = -0.0,-0.65,-1.3,-1.65,-2.00 and -2.35. 
The multi-component gas code SPEDQ was used for all 
gaseous models in the present paper. The power law IMFs 
were approximated by dividing the complete mass range into 
equal logarithmic steps. Alternative methods of discretiza- 
tion were also tried for certain cases in order to confirm 
the validity of the results, such as approximating the IMF 
using equal total masses in each of the components. The 
ranges of stellar masses {m max ,m m in) used in this paper 
are (1.0, 0.1), (2.0, 0.1) and (3.0, 0.1). The reason why higher 
values of m max were not used is that it is customary to sup- 
pose that a cluster would be largely depleted in heavier stars 
by the time gravothermal oscillations manifest ( |Kim, Lee fc| 
Goodman|19 98). In Sec|4j however, we briefly discuss a pos- 
sible exception. 

The critical value of N (N cr it) at which oscillations in 
the central density (p c ) first appeared (as N increased) was 
determined (correct to ten percent). The obtained values of 
N cr it in units of 10 4 are given in Table [I] 



sider the values of Mtot at N cr it {Merit)- The values of M cr it, 
for the models considered in the present paper, are given in 
Table[2] The values of M cr it in Tableware approximately the 
same for fixed mmax and the values vary much less with a 
than N cr it. Thus the conclusion of|Murphy ct al ( 1990 j) , who 
considered only evolved IMFs with fixed m m ax, appears also 
to apply to pure power law IMFs with fixed m ma x. What 
has been added in the study in the present paper is that the 
value of Merit also has a strong dependence on m max . 

We can compare the dependence on m ma x in Table [2 
with the results of the two-component models of |Breen fc 
Heggie (2012) if we fix the stellar mass of the light compo- 
nent in that earlier paper. This is done in Appendix |A| (see 
Table A2 I . In Table A2 there is a clear trend of increasing 
Merit with increased stellar mass of the heavy component for 
fixed total mass ratio (i.e. moving up through one column of 
Table A2|. Therefore the trend of increasing M cr u with in- 



creasing stellar mass range (or stellar mass ratio), for fixed 
nimin, seems to be a common feature of multi-component 
systems. It is also worth noting that for two-component sys- 
tems with fixed stellar mass ratio (see Appendix [A] Table 
A2 where we consider values along a given row) there is a 
trend of increasing Merit with decreasing total mass ratio. 
As decreasing total mass ratio for two-component systems is 
the analogue of decreasing a, these systems have the same 
stability trends as the multi-component systems in T able [2] 
Now we will attempt to extend the interpretation of [Breen] 



& Heggie (20121 from two-component to multi-component 



models in a way which also accounts for the dependence on 



1 SPEDI is a multi-component gas code which was ini- 
tially based on a formulation by |Louis fc Spurzem| (|1991[) and 
was subsequently further developed by |Spurzem fc Takahashi| 
(|1995|. Further information regarding SPEDI is available at 
http: / / www. ari . uni-heidelberg.de /gaseous-model / 



Breen & Heggie (20121 first argued that for the two- 



component case, the dynamics of the system were dom- 
inated by the heavier component. The reasoning behind 
this emphasis on Ns is that the heavier component con- 
centrates within the central region where it behaves like 
a one-component system deep in the potential well of the 
low mass stars, and it can exhibit gravothermal instability 
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like the central part of a one-component systerrj^] |Breen fc 
Heggie ( 2012[ ) showed that N2, the number of heavy stars, 
does indeed provide an approximate criterion for the onset 
of instability, and found the critical value of N2 to be of or- 
der 2000. For the case of multi-component models, however, 
it is unclear if and how the system can be divided into a 
heavy and a light component. Nonetheless, it may still be 
expected that the heavier stars may be more important to 
the dynamics of the system and to the onset of gravothermal 
oscillations. 



Table 3. Critical value of N e r in units of 10 4 . 



Breen & Heggie (20121 then gave an alternative mea- 
sure of the importance of the heavy component, in terms of 
what they called the "effective" number of heavy stars, and 
this is a concept that is more readily adapted to the case 



of several components. Breen & Heggie (20121 argued that, 



as the light component acts as a kind of container for the 
heavy component, it was the overall mass of this container 
(i.e. the total mass in the light component) that was the 
important factor and not the stellar mass of the light com- 
ponent. Therefore it could be replaced by an equal mass of 
stars of the massive component, giving rise to the idea of the 
effective number of stars N e f. This was defined as follows 



A, 



Mi + M 2 



(1) 



They also defined a modified half-mass relaxation time scale 
{tef,rh) by using N e f in place of N in the standard formula 
for the half-mass relaxation time. They used this effective 
relaxation time to modify and improve upon a stability cri- 



terion suggested by Goodman ( 19931; see Section 2.3 of the 
present paper. It is worth pointing out that N e f itself can 
be used as an approximate stability condition for the two- 
component models of Breen & Heggie (2012) (sec Appendix 
[A] ). It is this form of condition for the stability of two- 
component systems which we shall attempt to generalise to 
multi-component systems. 

Again it is not immediately clear what the most appro- 
priate extension of N e f to the multi-component case should 
be, as the appropriate definition of heavy stars is not clear. 
However, the result is much less sensitive to our choice than 
the number of heavy stars, A2. One simple approximate way 
to define a heavy star in this context is simply any star with 
a stellar mass of ~ m max . This would lead to the definition of 
N e f as m Itot for the multi-component model and no change 
in the definition for a two-component model (equation |l])). 
This value is given in Table|3j for the multi-component mod- 
els considered in this paper. We find that there is much less 
variation among the critical values of N e f than for M cr it, 
especially for varying m max . Nevertheless, the same trend 
of increasing M cr u with decreasing a as seen in Table [2] is 
still present in Table [3] in the form of increasing N e f. 

Now we will discuss a possible explanation for the in- 
crease in the critical value of N e f with decreasing a in Table 



2 Note that, in this picture, the gravothermal instability of the 
system is essentially confined to the massive stars; it is in the 
centrally concentrated massive stars that the temperature inver- 
sions occur which drive the expansion phase of the gravothermal 
oscillations. Around TV = N cr u, the light stars always exhibit a 
normal temperature gradient, and if their heat capacity is nega- 
tive this does not lead to instability. 
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[3] The idea behind N e f is that a multi-component system be- 
haves in approximately the same way as a single-component 
system consisting of N e f stars of stellar mass m max . As 
systems with higher a have more stars with stellar mass 
~ m max than systems with lower a, we may expect that 
the approximation with the one-component system is better 
for higher a than for lower a. Therefore, it is not surprising 
that in Table [3] it is the systems with a = that have the 
closest critical values of N e f to the critical value of N for a 
one-component system, i.e. about 7000. 

We can take our discussion further by considering sys- 
tems with fixed N e f. For systems with fixed N e f, as a de- 
creases there is an increasing number of light stars (stars 
with stellar mass not « m max ). As the number of light stars 
increases and the number of heavy stars decreases the two- 
body relaxation time increases, as it becomes increasingly 
dominated by the light component. According to Henon's 
principle ( Henon|[l975 1 the rate of energy generation in the 
system is regulated by two-body relaxation, and therefore 
there is a lower rate of energy generation as a decreases. 
Regardless of the value of a, the average mass in the core is 
approximately m max , and the lower rate of energy genera- 
tion can be met by a core of lower density. Thus for lower 
alpha there is a smaller density contrast between the core 
and the mean density in the system. This would imply that 
the stability would increase with decreasing a, as we indeed 



2.3 Goodman stability parameter 

Goodman ( 1993 I suggested the use of the quantity 



(2) 



^ _ Etot/t r h 
E c I ire 

as a stability indicator, where log 10 e ~ — 2 is the stability 
limit below which the cluster becomes unstable. Here E to t is 
the total energy, E c is the energy of the core, t rc is the core 
relaxation time and t r h is the half mass relaxation time. A 
condition of this type has been supported for two-component 



models by Kim, Lee & Goodman ( 1998 ) who studied Spitzer 



stable models using a Fokker-Planck code and by |Breen fc| 
Heggie (20121 who studied Spitzer unstable models using a 
gas code. However, |Breen fc Heggie (20121 also introduced 



a modified definition of e (see below) because the definition 
given in equation |2| was found to yield a critical value which 
varied with total mass ratio and stellar mass ratio. Also, the 
critical value at which the instability appears is somewhat 
different in|Breen & Heggie (2 012[ ) from that in Kim, Lee & 



Goodman (1998). 



For the multi-component models studied in the present 
paper, the values of log 10 e are given in Table [4] The values 
of log 10 e (based on the original definition, i.e. equation |2|) 
range from —2.26 to —2.56 and there is a decreasing trend 
with decreasing a. 
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Table 4. Critical value of 1 



Table 5. Critical value of J 



?10 £ 2 



(3,0.1) 
(2,0.1) 
(1,0.1) 







-0.65 



-1.3 



-1.65 



-2.00 



-2.35 



(m„ 



iin )\f 



-2.26 
-2.26 
-2.26 



-2.38 
-2.32 
-2.30 



-2.54 
-2.48 
-2.39 



-2.61 
-2.53 
-2.40 



-2.57 
-2.58 
-2.43 



In equation Q, 
0.34ag 

~ G 2 



t rc and t r h are defined by 



m c p, 



In A 



(-2.62) 

-2.56 

-2.42 



(3) 



and 



0.138iV5r^ 



(4) 



(Gm) 5 In A 

where m c and a c are the mass density weighted averages 
over all the components in the core. However, the definition 
of t r h does not take into account the mass spectrum. |Breen| 
& Heggie ( 2012 ) have shown that for two-component models, 
modifying the definition of t r h to take into account the mass 
spectrum leads to an improved stability condition. As has 
already been shown in Section [2. 1| we can construct a value 
N e f which provides an approximate stability condition. Now 
we will use this value to define a new half mass relaxation 
time defined as 



0.138^^ 



ef,rh 



(5) 



(Gm m „)2 In A 

We use this in place of t r h in equation ([2]) to de fine the new 
stability parameter £2 as in |Breen fc Heggie| (|2012), The 
values of this parameter are given in Table |5| The variation 
of log 10 €2 in Table [5] is of comparable magnitude to the 
variation of log 10 e in Table [4] but the values in Table [5] 
are more consistent with that of a one-component model 
(log 10 e = log 10 £2 = —2) than those in Table [4] 

For the one-component model the definitions of e and 
£2 are identical. Also, it is worth noting that, for two- 
component systems with extremely small amounts of the 
heavy component (relative to Mtot), tef.rh as defined in 
equation (|5| is not a suitable approximation to the relax- 
ation time. This is the case for the models considered by 
Kim, Lee & Goodman (19981, and so the extension to €2 



would not have been necessary or useful in the context of 
their paper. 

The logarithm of the Goodman stability parameter (or 
our somewhat more consistent modified version) seems to 
have a particular value approximately —2 at the stability 



boundary (Kim, Lee & Goodman (19981, Breen & Heggie 



(20121 and the present paper). However, it is not known if 
£ (or £2) can be predicted for a particular IMF without car- 
rying out numerical simulations, which limits its usefulness. 
In contrast N e f has the advantage that it can be easily cal- 
culated before carrying out numerical simulations, and also 
provides an approximate indication of the stability bound- 
ary. 



3 DIRECT N-BODY SIMULATIONS 

In order to validate the results obtained from the gas code 
SPEDI a series of TV-body runs were carried out. The direct 
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TV-body simulations in the present paper were conducted 
using the NBODY6 code ( jAarseth|20"03"| |Nitadori fc Aarseth] 
2012[ ). As the IMF's with a = in Table |1] have the lowest 
values of N cr it (2.0 x 10 4 ), for this value of a TV-body runs 
were carried out for the mass ranges (1,0.1) and (2,0.1) . 
The case with parameters a = —1.3 and (m ma x, fnmin) = 
(1.0,0.1) was also chosen because it has a higher value of 
N cr u (3.5 x 10 4 ) than for a — 0, although the value is still low 
enough to make it suitable for direct A-body simulations. 

For the case of a = 0.0 the values of N cr u are the same 
(see Table [T| regardless of the stellar mass range. For the 
two mass ranges chosen, there were no signs of gravothermal 
behaviour in the A-body runs with TV = 8fc or TV = 16k. 
The first clear sign of gravothermal behaviour occurs with 
N — 32k for both chosen mass ranges. For the stellar mass 
range (1.0, 0.1), the three panels of Fig.[l]show, respectively, 
(i) the evolution of the core radius r c , (ii) an example of a 
single cycle of gravothermal oscillation in the post-collapse 
evolution, and (iii) evidence of the gravothermal nature of 
the oscillation for the 32fc run. The same graphs for the 32fc 
run with stellar mass range (2.0,0.1) are given in Fig. [2] 

For the case a = —1.3 and stellar mass range (1.0,0.1) 
the value of Ncrit is 3.5 x 10 4 (see Table [l]). No gravother- 
mal behaviour was seen in the TV-body runs with TV = 8k, 
12k and 32fe for this set of conditions. The first signs of 
gravothermal behaviour occurred in the 64A: run as would 
be expected from the above value of N cr it obtained from 
SPEDI. The same three graphs shown for both the a — 0.0 
cases (see previous paragraph) are plotted for the 64fc run 
inFig|3] 

In the graphs of r c for all cases (see Fig. [l]top, Fig. [2] 
top and Fig. [3] top), behaviour can be seen which is qual- 
itatively similar to gravothermal oscillations (see Makino| 
( p96l),|Takahashi fc Inagaki] ( fl995| > , |Heggie fc Giersz | ( |2009f 
and Breen & Heggie (2012 1). For the case of a = 0.0 with 
the mass range (1.0,0.1) (Fig. [l]top) one oscillation can be 
seen between 3830 and 4570, and another between 5970 and 
6560. During each of these oscillations r c changes by more 
than a factor of 10. Similarly for the case of a = 0.0 with the 
mass range (2.0,0.1) (Fig. |5]top) an oscillation in r c can be 
seen between 6800 and 7950, and part of an oscillation can 
also be observed after 8610. In this model the change in r c 
is about a factor of 10. Finally for the case of a = —1.3 with 
the mass range (1.0, 0.1) (Fig.[3]top) an oscillation in r c can 
be seen between 4800 and 5600, and part of an oscillation 
can also be observed after 7400. The change in r c is about 
a factor of 10, which is similar to the change in r c for the 
a — —1.3 runs. 

Now we consider the physical nature of these oscilla- 
tions, which could in principle be driven by sustained binary 
activity or by gravothermal behaviour. A sign of gravother- 
mal behaviour is that the binding energy of the binaries 



remains roughly constant during times of expansion (McMil- 
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|lan fe Engle||1996| ). This is because the expansion phase of 
a gravothermal oscillation should be driven by the core ab- 
sorbing heat from the rest of the cluster rather than by en- 
ergy generation. At core bounce, where p c reaches a local 
maximum, there is an increase in binary activity, and enough 
energy is produced to halt and reverse the collapse. This be- 
haviour is particularly clear in Fig. [3] (middle) where there is 
an initial increase in relative binding energy of the binaries 
coinciding with core bounce and the initial expansion. A bi- 
nary escapes, and then there is a period of expansion during 
which the relative binding energy of the binaries remains 
nearly constant (from 4850 to 5000). Towards the end of the 
oscillation there is renewed binary activity corresponding to 
the next core bounce. There is also binary activity at other 
times during the oscillation, but it has no discernible effect 
on the evolution of r c . Fig. [I] (middle) is also a good example 
of gravothermal behaviour. Mild binary activity continues 
after core bounce (from t = 3860 to t = 3930) but expan- 
sion continues thereafter for a period. However in Fig. [2] 
evidence of gravothermal behaviour is more ambiguous. We 
will discuss the case of Fig. [2] in detail in the last paragraph 
of this section. 

The cycles of p c vs the core velocity dispersion v\, as 
seen in Fig. [I] (bottom), Fig. [5] (bottom) and Fig. [3] (bottom), 
are believed to be a sign of gravothermal behaviour ( |Makino] 
1996). During these cycles, the temperature is lower during 



the expansion where heat is absorbed and higher during the 
collapse where heat is released. The velocity dispersion in 
Fig [2] (bottom) and Fig. [3] (bottom) has been smoothed to 
make the cycle clearer. These cycles are similar to the cycles 
found by|Makino|(|1996[ ) for one-component models and by 
Breen & Heggie 12012 1 for two-component models. 



The gravothermal nature of the behaviour is clearer in 
Fig. [l] (where {m max , m mtn ) = (1.0,0.1)) than in Fig. [2] 
(where {m max , m m in) = (2.0,0.1)). This is perhaps surpris- 
ing as both cases have the same value of N cr u as found with 
SPEDI. We will now discuss a number of possible reasons for 
this apparent difference in behaviour. Firstly, as the values 
of N cr u are only correct to 10%, it is possible that in reality 
the values could differ by up to 4 x 10 3 . Secondly, another 
issue is that in the gas model the mass function is discre- 
tised, resulting in a difference between the mass of the heavy 
component (mio) and m max . mio is about 14% percent less 
than m ma x for the stellar mass range (2.0, 0.1) and 11% for 
(1.0,0.1). ft is argued in Appendix [A] that the stability of 
a system will increase if the average stellar mass inside the 
core is increased (while keeping the stellar masses outside 
the core approximately the same). This would imply that in 
both cases the values of N cr it found with the 10-component 
models are underestimates for the onset of instability, and 
that the true value of N cr it for (2.0,0.1) is slightly higher 
than for (1.0,0.1). Finally, the fact that the gravothermal 
nature of the behaviour is clearer in one run might simply 
be a stochastic effect. 

For the purposes of this paper, we have not considered 
the evolution of multi-component systems in the regime N < 
N cr it. In Spitzer-unstable cases, there is no reason to doubt 
that this is characterised by mass-segregation, followed by 
post-collapse expansion powered by binary evolution as in 
the much smaller N-body models considered long ago by|van| 



4 SUMMARY AND DISCUSSION 

The focus of this paper has been on the conditions for the 
onset of gravothermal oscillations in multi-component sys- 
tems. We have investigated power law IMFs with different 
exponents and three different stellar mass ranges (3.0,0.1), 
(2.0,0.1) and (1.0,0.1). A multi-component gas code has 
been used to obtain the values of N cr u- In order to verify 
the validity of the results direct TV-body runs were carried 
out on appropriately chosen cases. The values of N cr it found 
ranged from 2 x 10 4 to 10 s , and varied with a and the stellar 
mass range. 

Motivated by |Murphy et al| ( [1990] ), who found that the 
total mass of the systems they studied could be used as an 
approximate stability condition, the value of Merit (the total 
mass of the system at N cr u) f° r each system was calculated 
(see Table [2). While for a fixed mass range Merit does pro- 
vide an approximate stability condition, the value of M cr it 
varied by roughly a factor m max . Merit can also be used as 
an approximate stability condition for the two-component 
models of |Breen fc Heggie ( 2012| so long as the stellar mass 
ratio is fixed (see Appendix |A[). 

In order to find a more general stability condition we 
applied an extension of an idea first employed in |Breen fc] 
Heggie ( 2012 \. They used a quantity called the effective par- 
ticle number (N e f)- The value was useful because the two- 
component system that was being considered was expected 
to behave in roughly the same manner as a one-component 
system with N e f stars. In the present paper this idea has 
been extended to multi-component systems. The values of 
N e f for the multi-component models in this paper are given 
in Table [3] The variation in Table[3]is significantly less then 
that in either Table [l] or Table [2j A stability condition of 
N e f ~ 10 4 covers most of the values of Table |3| and indeed 
the two-component models of |Breen fc Heggie ( 2012J (see 
Table |A3l. 



The Goodman Stability Parameter was also tested for 
the multi-component case (see Table |4j). The critical val- 
ues m Table |3 were found to be lower than the value for a 
one-component model (log 10 e = —2) and also varied with 
a and, to a much lesser extent, with m max . By modifying 
the Goodman Stability Parameter using a slightly different 
definition for the half-mass relaxation time (based on the ef- 
fective particle number) a critical value was found which was 
more consistent with the critical value for a one-component 
model (see Table |5b. 

Goodman] ( |1987[ ) used a gas model to find the value of 



N cr u (= 7000) for a single component system. Technically 
what he showed was that steady post-collapse expansion was 
possible in a gas model for all N, but that it was unstable for 
N > 7000. While the gas model used by|Goodman|(|1987|) is 



similar in form to the model used here and in|Breen fc Heggie 



(2012). there are two notable differences. Firstly 



Goodman 



( 1987) used a larger energy generation rate than the one used 
here. Secondly, the parameter of the coulomb logarithm that 
was used was A = 0.4. A value of A = 0.4 ( |Spitzer|1987) was 
a reasonable choice at the time, but it has since been shown 
that A = 0.11 is a better choice for a single-component model 
(Giersz fc Heggie||1994|. (For multi-component models the 



Albada (19671, Aarseth (19681 and many more since 



value of A = 0.02 was found to provide a good fit (Giersz & 
Heggie ||l996 )). These two differences affect the stability in 



opposite ways: by arguments similar to those given in the 
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last two paragraphs in Appendix [X] a larger energy gener- 
ation rate will increase stability, whereas a larger value of 
A tends to reduce stability. For example for N = 7000 with 
A = 0.4 the increase in the relaxation rate is 20% compared 
with A = 0.11. 

In the present paper, we have made the assumption that 
multi-component systems will be depleted in stars with stel- 
lar mass greater than 3Mq. This neglects the possibility of 
systems containing a population of stellar mass Black Holes, 
which would require a value of m max about an order of mag- 
nitude greater than what is considered here. These systems 
are outside the parameter space studied by |Breen fc H eggie 
(20121 and Kim, Lee & Goodman (19981, as the total mass 



ratio is lower than the range considered by |Breen fc H eggie 
(20121 and the stellar mass ratio is higher than the values 



considered by Kim, Lee & Goodman (19981. The onset of 



gravothermal oscillations and the more general evolution of 
systems containing a population of stellar mass black holes 
are the topics of the next paper in this series. 

To conclude, a stability condition of N e f ~ 10 4 does 
apply to the multi-component systems in this paper and 



the two-component systems of Breen & Heggie (2012). This 



condition is expected to apply to any multi-component sys- 
tem provided that there is a sufficient number of stars with 
stellar mass ~ m max . 
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APPENDIX A: THE TWO-COMPONENT CASE 
REVISITED 

The purpose of this appendix is to reconsider the results 



of Breen fc Heggie] ( |2012[ ) in terms of the effective par- 
ticle number N e f defined in equation [l] |Breen fc H eggie 



(2012) investigated gravothermal oscillation in a range of 



two-component models, specified by the stellar mass ratio 
^ and total mass ratio where m,2 (mi) is the stellar 
mass of the heavy (light) component and M2 (Mi) is the 
total mass of the heavy (light) component. For reference the 
values of N cr it for these models are given in Table [AT] which 



is similar to Table 2 in Breen & Heggie (2012 1. The differ- 



ence is that the data have been rearranged to compare more 
closely to the ar rang ement in the present paper. Thus the 
columns in Table A 1 are arranged in order of decreasing ^ 
as this is the analog for two components of decreasing a. 

Following the approach in the main part of this paper, 
we will firstly consider the values of M cr it- In order to con- 
sider M cr u we need to specify the mass unit, and to make 
the results comparable with the multi-component models in 
the main part of the present paper mi has been fixed at 
O.IMq. The values of M C rU for the two-component models 
in Table [AT] are given in Table |A"2| For comparison the value 
of M cr it for a one-component model would be 0.07 x 10 4 
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Table Al. Critical value of N (N cnt ) in units of 10 4 



m 2 \ M2 
m-i * Mi 


1.0 


0.5 


0.4 


0.3 


0.2 


0.1 


50 


18 


30 


33 


42 


55 


100 


20 


8.5 


13 


15 


18 


22 


36 


10 


5.0 


7.2 


8.2 


10 


12 


22 


5 


2.8 


1.0 


1.6 


5.4 


7.0 


10 


4 


2.4 


3.5 


3.8 


4.6 


5.5 


8.5 


3 


2.0 


2.8 


3.2 


3.6 


1.1 


6.0 


2 


1.7 


2.2 


2.3 


2.6 


3.0 


3.8 



Table A2. Critical value of M crit in units of lO^Mp 



©■ 



The value 
■f is 0.07 x 



of mi is fixed at O.1M0. For reference the value of M, 
1O 4 M0 for-^jj- = 0, which is obtained from the result of|Goodman 
1 1987 1 for a one-component cluster with mi = O.IMq 



m 2 \ M 2 
mi ' Mi 


1.0 


0.5 


0.4 


0.3 


0.2 


0.1 


50 


3.529 


4.455 


4.583 


5.427 


6.574 


10.978 


20 


1.619 


1.902 


2.059 


2.305 


2.614 


3.940 


10 


0.909 


1.029 


1.104 


1.262 


1.412 


2.396 


5 


0.467 


0.545 


0.596 


0.662 


0.808 


1.078 


4 


0.384 


0.467 


0.484 


0.556 


0.629 


0.912 


3 


0.300 


0.360 


0.395 


0.425 


0.495 


0.639 


2 


0.227 


0.264 


0.268 


0.294 


0.327 


0.398 



(using m = O.1M0 ). Th is is significantly lower than any of 
the values in Table 



A: 



Me 



For fixed values of — the value of 

mi 

lcrit varies by factors of up to « 3 between = 1 and 0.1. 
Therefore for fixed ^ > M crit does provide a rough stability 
condition. However for fixed j^-, the variation in M cr u is 
a factor of ~ 15 — 30. The variation of M cr n with varying 
rri2 resembles the variation of M cr u with varying m max in 
Tabled 

We will now consider the values of N e f for the two- 
component systems; these are given in Table |A~3| For com- 
parison the critical value of N e f for a one-component model 
is the same as its value of N crit , which is 0.7 x 10 4 . The 
values in Table |A3| vary much less then those of M cr i t in 
Table [A2| although, as pointed out in Section [2. 1| N e f can 
be interpreted as a measure of the total mass of the system 
in units of m.2. A stability condition of N e f ~ 10 4 or slightly 
more covers, within a factor 2 at most, the values of Table 
]and indeed Table [3] 



A2 



3 In Table 

Given that M cr u 
the values in Table 



for fixed i Merit increases with decreasing 
for 

"4 1 



observed in Table 



,,, 1( - — is significantly lower than any of 
|A1| one may wonder if that increasing trend 



Al 



continues below 



0.1. This is a topic 



considered in detaiTm the next paper of this series. 



Table A3. Critical value of N e f in units of 10 4 



m.2 \ M 2 
mi \ Mi 


1.0 


0.5 


0.4 


0.3 


0.2 


0.1 


50 


0.71 


0.89 


0.91 


1.09 


1.31 


2.20 


20 


0.81 


0.95 


1.02 


1.15 


1.31 


1.97 


10 


0.91 


1.03 


1.10 


1.26 


1.42 


2.40 


5 


0.93 


1.09 


1.19 


1.32 


1.62 


2.16 


4 


0.96 


1.17 


1.21 


1.39 


1.57 


2.28 


3 


1.00 


1.20 


1.32 


1.42 


1.65 


2.13 


2 


1.13 


1.32 


1.34 


1.47 


1.64 


1.99 



All of the trends in Table |A3| may be understood if 
we consider the reasoning behind the use of N e f as an ap- 
proximate stability condition. The basic idea is that the 
multi-component system in question evolves in a similar way 
to a one-component system of N e f stars with stellar mass 
77i2 ■ This requires that the half mass relaxation timescale of 
the multi-component system is similar to that of the one- 
component system with which we are comparing it. We as- 
sume this to be true if the heavy component amounts to 
a significant faction of the total mass within the half-mass 
radius rn, which is certainly not the case as tends to 0, 
i.e. on the extreme right of Table |A"3| 

We now consider with more care how the two- 
component system actually differs from the corresponding 

and 2a ar e 

mi 

varied. For fixed 



one-component system as the parameters ^ 

. decreases the relaxation process 

is increasingly dominated by the light stars. This leads to 
the system behaving more like a one-component system of 
M ^ stars as opposed to a one-component system of N e f 
stars, and this increases the half-mass relaxation time. As 
the rate of two body relaxation becomes slower the core be- 
comes larger (relative to Th', see the discussion of Henon's 
Principle in Section 2.2 1 as it can produce the required en- 



ergy at a lower mass density (as the average stellar mass in 
the core remains the same, roughly m-£). Because gravother- 
mal instability depends on a high density contrast within the 
system, it would be expected that stabil ity w ould increase 
decreases, as can be seen in Table 



M-, 
M L 



A3 



Now let us consider the case of fixed If we consider 

Mi 

the post collapse evolution of series of systems with fixed M2 
and m2, as ^ decreases the tendency towards mass segre- 
gation becomes weaker. Therefore the half mass radius of the 
heavy component (rhfi) is smaller compared to for larger 

— than for smaller — . It follows that the mass density 
of the heavy component (within ^,2) is smaller for smaller 

— than for larger — . The relaxation time of the heavy 

mi mi ^ 

component within its half-mass radius ru.i itrh,2) decreases 
with increasing mass density. This leads to the conclusion 
that the relaxation time within the heavy component in- 
creases with decreasing The energy flux in the heavy 
component, which we are assuming regulates the rate of en- 
ergy generation, is of order (where E2 is the energy of 
the heavy system). Therefore, as ^ decreases so does the 
energy flux, which results in a lower rate of energy genera- 
tion. The lower rate of energy generation leads to a larger 
core (relative to rh) as the core can produce the required 
energy at a lower mass density. Thus it would be expected 
that stability (as measured by N e f) would increase as ^ 
decreases, and this is what is observed in Table I A3I for most 



values of ^t 2 - • However, the trend of increasing stability with 



decreasing ^ seems to disappear for small . Reasc 
this will be discussed in the next paper of this series. 



This paper has been typeset from a TjrjX/ P/TjrjX file prepared 
by the author. 



© 2010 RAS, MNRAS 000,[T}{7] 



8 P. G. Breen and D. C. Heggie 




1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 2 3 4 5 6 7 8 9 

log p log p 

Figure 1. ./V-body run of a multi-component model with N = Figure 2. As Fig. [l] with N = 32k and a = 0.0, but 

32k, a = 0.0 and ( 

Trijriax , Tfimiri) — (1.0,0.1). Top: log v c vs time {wimax , ^min) — (2.0,0.1). The velocity dispersion has been 

(N body units) over the entire run. Middle: logr c vs time (TV smoothed to make the cycle clearer, 

body units) over-plotted with the relative binding energy of bi- 
naries. The plot is over the period of a gravothermal oscillation 
which occurs between 3830 and 4570. Bottom: p c vs showing 
the gravothermal nature of the cycle over the same time period 
as the middle plot. 
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-2.6 - 
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log v 2 vs log p over time 4750 to 5600 
| 1 1 1 1 




-0.25 1 ' ' 1 1 1 ' 1 

1.6 1.8 2 2.2 2.4 2.6 2.8 3 

log p 

Figure 3. As Fig.[l] with (mmax , rn min ) = (1.0,0.1), but N = 
6Ak and a = —1.3. The velocity dispersion has been smoothed to 
make the cycle clearer. 
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